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ABSTRACT 


In this thesis, the impacts of transforming the coordinate system of an existing 
x-z mesoscale model to are analyzed and discussed as they were observed in 
three test cases. The three test cases analyzed are: A rising thermal bubble, a 
linear hydrostatic mountain, and a linear nonhydrostatic mountain. The methods 
are outlined for the transformation of two sets (set 1, the non-conservative form 
using Exner pressure, momentum, and potential temperature; and set 2, the non¬ 
conservative form using density, momentum, and potential temperature) of the x-z 
Navier-Stokes equations to x-a^ and their spatial (Continuous Galerkin) and temporal 
(Runge-Kutta 35) discretization methods are shown in detail. For all three test cases 
evaluated, the x-a^ models performed worse than their x-z counterparts, yielding 
higher RMS errors, which were observed predominantly in intensity values and not 
in placement of steady state features. Since the models did converge to a fairly 
representative steady-state solution the results found by this project are promising, 
even though they did indicate that x-Cz coordinates are not as accurate or efficient 
as x-z coordinates. With further fine-tuning of the model environment, these issues 
could be made minimal enough to warrant their utility with semi-implicit methods. 
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I. 


INTRODUCTION 


Numerical Weather Prediction (NWP) models are the work horse of mod¬ 
ern atmospheric condition forecasting and as such have been the focus of numerous 
studies, aiming to develop more accurate models with longer deterministic forecast 
periods. There are countless areas in which to focus this research (spatial and tem¬ 
poral discretization methods, non-reflecting boundary conditions, data assimilation, 
prognostic equations, physical parameterizations, etc....), and for this project the fo¬ 
cus will be on the transformation of the prognostic or governing equations from x-z 
coordinates to terrain following x-a^ coordinates using a specihc spatial (continu¬ 
ous Galerkin) and temporal (Runge-Kutta 35) discretization method. This thesis is 
necessary in order for future research to be able to evaluate and compare various 
coordinates systems while varying the temporal discretization methods, allowing for 
larger time steps while maintaining stability (i.e., semi-implicit methods). 

Various mature mesoscale models such the Coupled Ocean/Atmosphere Meso- 
scale Prediction System (COAMPS) [2] and the Weather Research and Forecasting 
(WRF) modeling system [3] use a variation of the x-a coordinate transformation. By 
studying the works of Gal-Chen and Somerville [4] and analyzing the transforma¬ 
tion methods used in COAMPS and WRF, a similar x-a^ coordinate transformation 
was applied to the governing equations of interest employing continuous Galerkin 
techniques. 

The governing equations selected for this thesis consist of two forms of the 
Navier-Stokes equations. The Navier-Stokes equations, along with their variations, 
form the most widely used and accepted sets of equations for numerically resolving 
atmospheric flow. The first specihc formulations of the equation sets selected (set 
1) was the non-conservative form using Exner pressure, momentum, and potential 
temperature, which is used in the operational NWP model COAMPS [2]. The second 
formulation chosen (set 2) was the non-conservative form using density, momentum. 
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and potential temperature. The operational NWP model WRF [3] also uses set 2, 
but in a conservative form. Building on the work of Giraldo [1] who by implementing 
continuous Galerkin techniques, developed a 2-D (x-z slice) mesoscale model using 
Non-Hydrostatic Equations (Euler and Navier-Stokes Equations), the original con¬ 
struct was transformed from x-z coordinates to x-a^ coordinates to test the impacts 
on resolving atmospheric motion in a continuous Galerkin (GG) framework. Gur- 
rently, most operational non-hydrostatic models use hnite difference (ED) methods 
(i.e., structured grids), which then rely on x-a in order to resolve atmospheric flow 
in the presence of terrain. GG methods, on the other hand, can use various types 
of grids (unstructured grids, x-z grids, x-a^ grids, etc....). Eor this reason, the GG 
method is well positioned to judge the effects of the coordinate system on the solution 
accuracy and efficiency, which is the goal of this thesis. 

After the modihcations were made to the model, three test cases were run: ris¬ 
ing thermal bubble, linear hydrostatic mountain, and linear non-hydrostatic moun¬ 
tain. The numerical solutions were either evaluated against other model solutions 
(case 1) or the analytic approximations (case 2 and case 3) using root mean squared 
error and normalized momentum flux. The resultant data was also compared to 
the unmodihed solutions. Additionally, the initial conditions for the test cases were 
pre-dehned for each case, maintaining uniform initial conditions from which both 
coordinate systems numerical solutions can be compared. 

With x-(Tz coordinates that are proven to function properly using fully explicit 
time integration, future research will be able to evaluate x-a^ coordinates using semi- 
implicit methods. NWP models are already taking advantage of semi-implicit time 
integration methods, which optimize the horizontal and vertical resolutions and their 
associated time sets while maintaining stability. The 2-D semi-implicit method (x- 
z formulation) can use very large time-steps sizes since the Gourant-Priedrichs-Lewy 
(GEL) condition is no longer constrained by acoustic and gravity waves; the penalty is 
that a global 2-D implicit problem must be solved. In contrast, the 1-D semi-implicit 
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x-az formulation can use large time-steps, but it must adhere to the CFL condition 
due to gravity and acoustic waves in the horizontal; the advantage is that only a 1-D 
matrix problem needs to be solved. A model using x-az can exploit semi-implicit 
methods along the vertical (a) direction, which is the inspiration of this thesis topic. 
Both COAMPS and WRF are currently using semi-implicit methods (only along a), 
but are constrained by FD spatial discretization methods. 

Though only explicit time integration was used in this thesis, it is now possible 
to observe if the implementation of x-az coordinates signihcantly improves or dimin¬ 
ishes the solution of the Navier-Stokes equations over x-z coordinates when using a 
CG framework. Giraldo [5] has already developed semi-implicit methods for the 2-D 
model in the x-z framework and tested their impacts. With the model coordinates 
transformed to x-az, future research will be able to extent the time integration to 
semi-implicit methods for x-az and compare the results with the semi-implicit results 
derived using x-z coordinates. The relevant governing equations for this project, the 
coordinate transformation theory, and the discretization methods are discussed in 
Chapter II. The application of the coordinate transforms are discussed in Chapter 
III. The three test cases are explained and outlined in Chapter IV. A discussion and 
interpretation of the resulting impacts from the transformation is in Chapter V. The 
conclusions found and recommendations are presented in Chapter VI. 
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II. 


BACKGROUND 


A. GOVERNING EQUATIONS 


1. Equation Set 1: Navier-Stokes Equations with 
Exner Pressure 

The first set of equations chosen was the Navier-Stokes eqnations that uses 
Exner pressnre, of which there has been extensive amonnts of documented work (i.e. 
this is the formulation used in COAMPS) for comparison. This set can only be writ¬ 
ten in non-conservative form and consists of a system of three eqnations. The hrst 
equation is the pressure tendency eqnation: 


Ott R ^ ^ 

——h u ■ Vtt H-ttV ■ u = 0 

at c„ 


( 2 . 1 ) 


71 = 



JL 

Cp 


where tt is Exner pressure, u represents the velocity held (u, w), R is the gas con¬ 
stant, Cv is the specihc heat for constant volume, Cp is the specihc heat for constant 
pressure, P is pressure, and Pq is pressnre at the snrface. The second equation is the 
momentum equation: 


du 

'm 


+ u ■ Vu + CpOVn 


—gk + 


( 2 . 2 ) 


where 9 is potential temperatnre, g is the gravitational constant, /c is a vector (0,1)^, 
and fi is the dynamic viscosity. The third eqnation is the thermodynamic energy 
eqnation: 
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+ M = /iV"^ 
at 


(2.3) 


For the scope of this project, only inviscid flow will be considered (i.e. /x = 0) and the 
equations further reduce to the Euler equations which will be used in the following 
sections. 


2. Equation Set 2: Navier-Stokes Equations with Den¬ 
sity 

The second set of equations chosen was a version of the Navier-Stokes equa¬ 
tions that is now used in contemporary NWP models (i.e. this is the formulation 
used by WRF) and uses density, momentum, and potential temperature as the pri¬ 
mary state variables. Unlike set 1, this set can be written in both conservative and 
non-conservative form. Though neither form can conserve energy, using the non¬ 
conservative form can still conserve mass and more sophisticated time integration 
strategies can be used [6]. Thus for this thesis the non-conservative form will be used, 
which consists of a system of three equations. The hrst equation is the mass equation; 

^ + V . (pu) = 0 (2.4) 


where p is density. The second equation is the momentum equation: 


— + u - Vu + -VP = —gk + uV^M 
at p 


P 


Po 


m 


(2.5) 
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where P is pressure, Pq is pressure at the surface, and 7 is —. The third equation is 
the thermodynamic energy equation: 

+ = ( 2 . 6 ) 


Similar to set 1, only inviscid flow will be considered (i.e. /i = 0) and the equations 
further reduce to the Euler equations which will be used in the following sections. 

B. X-Z TO X-dz COORDINATE SYSTEM TRANSFORM 

1. Gal-Chen and Somerville 

In 1975, Gal-Chen and Somerville took the anelastic approximation of the 
Navier-Stokes Equation (in the cartesian form) and transformed the coordinated sys¬ 
tem to sigma-z coordinates [4]. An expanded set of prognostic equations was used 
for Gal-Chen and Somerville’s derivation and only the hrst three equations and the 
resulting transform will be used for comparison in this project. The hrst equation 
was the continuity equation: 


{Pou^),j= 0 . 


(2.7) 


where po is density, u is the velocity components (u, v, w )^, and j is the derivative 
operator where: 


d 
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The second equation was the momentum equation; 



(poti') + (pi>u‘u’),j = +S‘^p'g + T‘-'y . 


( 2 . 8 ) 


where is the Kronecker delta, p' is the density perturbation, p' is the pressure 
perturbation, and is the eddy viscosity. The third equation used was the thermo¬ 
dynamic energy equation: 


St 


(poO') + (poO'u^),j = H^, 


(2.9) 


where 9' is the perturbation in potential temperature and is the eddy diffusion. 

Using equations 2.7 - 2.9, Gal-Chen and Somerville derived the following set 
of transformations [4]: 


X = x, y = y, 


_ _ H{z - Zs) 
{H-Zs) 


m _ dzsZ-H m _ dzsZ-H m _ H 

dx dx H — Zs' dy dy H — Zg' dz H — Zg 


' u ' 


V 

= 

, w , 



1 , 

0 , 


0 , 

1 , 


0 

0 




dzs z—H dzs z—H z—H 

y dx H—Zs ’ dy H—Zs ’ H—Zs ) 


V 

yw j 


where the variables with the ^ represent the variables that have been transformed, H 
is the height at the top of model space, and Zg is the height at the surface of the 
model. 





i) ii) iii) 

Figure 1. x-z to x-cr^ coordinate transformation: (i) traditional x-z coordinates, (ii) 
x-a^ coordinates, and (iii) x-cx^ coordinates mapped back to x-z space. 


These transformation functions are used to convert traditional x-z (see figure 
l.i) to x-az coordinates (see figure l.ii), which when mapped back to x-z space are 
terrain following (see hgure l.iii). The derived inverse transformations are: 


.z(H-Zs). , 

x = x, y = y, z = [ - — -J + 2 ;, 


u 


V 

= 

L W , 



0 , 

dzs z—H 


0 , 

1 , 


0 

0 




_ _ _ _ dzs z—H H—Zs 

\ dx H ^ dy H ^ H j 


V 

yw J 


2. Basic Transformation Machinery 

This section will outline the basic equations used to transform the x-z coordi¬ 
nates of the Navier-Stokes Equations to x-a^, which are similar to the transformations 
derived by Gal-Chen and Somerville in the previous section, but in two dimensions. 
The hrst concept used was the total differential: 


dx 


dx dx 

Tz-dx + —dz 
ox Oz 
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da. 


daz 

dx 


dx + 


daz 

dz 


dz 


where the'notation indicates the transformed variable. The two total derivatives can 
then be written as a system of equations: 


dx 
^ daz 


/ 

dx 

dx 

\ 


dx 

dz 



daz 

daz 


V 

dx 

dz 

) 


dx 
^ dz 


( 2 . 10 ) 


Using vector notation (^) the above system can be simplihed to: 


dx = Jdx 


( 2 . 11 ) 


where J is the Jacobian of the transformation. Considering that the velocity u is 
dehned by the transformation can then be manipulated to become: 


u = Ju 


( 2 . 12 ) 


Using basic linear algebra and assuming that J is invertible, which is true for an affine 
(one-to-one) mapping, the inverse transform for velocity can be written as: 


u = J 


(2.13) 


The next major mathematical concept derived for the transformation was the gradi¬ 
ent operator (V) which, as a system of equations, can be written as: 
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d d dx d daz 

dx dx dx daz dx 

d d dx d daz 

dz dx dz daz dz 


which when put into matrix form becomes: 


A ^ 


dx da;^ ^ 


( A ) 

dx 


dx dx 


dx 

A 


dx da^ 


d 

dz j 


\ dz dz j 


V j 


Similar to the total derivative, the gradient, using vector notation, was the defined as: 

V = J^V (2.14) 


The last major concept used for the coordinate transformation is the linear algebra 
identity: 


{ABf = 


(2.15) 


where A G and B G The matrices were multiplied together and then 

transposed and also rearranged, transposed, and then multiplied to show equality: 


{ABf = 


... Gt'l,r 


••• j y K,i ■■■ bn,i j 


T 


&i,i ••• hi^i 
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+ ••• + ••• ®1,1^1,/ + •••+ 

flm,1^1,1 + ••• + (lm,nbn,l ■■■ 0-m,lbl,l + ■■■ + Clm,nbn,l 


Ol, 1^1,1 + 

... + ai^nbn,! 

flm, 1^1,1 + • 

• • (^7n,nbn,l 

1 - 

+ 



Ctrn,nbn,l 



( h 

Ol,l 

b ) 

On,l 


/ 

Ctrn,! 



^ hi,i ... 

bn,l j 


y CH,n 

Clm,n J 

1 ^ 1,1 + •• 

• + 0-l,nbn,l 

®m,l^ 

'1,1 + 


+ •• 

■ + 0-l,nbn,l 


h,z + 

• “1“ 


3. Transformation Functions 

For the transformations from x-z coordinates to x-a^, the x and ^ coordinates 
from Gal-Chen and Somerville [4] were rewritten using new notation which will be 
seen through the remainder of the thesis: 


X = X, (Jz 


H{z - Zs) 

{H - Zs) 


(2.16) 


where x is now x and 2 : is Using Eqs. (2.16) their derivatives were evaluated, 
yielding: 


d(7z dzs az — H d(7z H 
dx dx H — Zs' dz H — Zs 


( 2 . 17 ) 
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which matched the derivatives found by Gal-Chen and Somerville [4], The derivatives 
were then applied to various forms of the Jacobian from Eq. (2.10), which yielded: 


J 


1 , 0 

dzs a,-H H 
\ dx H-Zs ’ H-Zs 


( 2 . 18 ) 





dzs ctz—H 
dx H — Za 



H 

H-Za 


(2.19) 


J-' 


/ 


1 , 


0 


I dZa Uz—H H — Za 

\ dx H ^ H 


( 2 . 20 ) 


{j-r 



dZa CTz—H 

dx H 


H-Za 

H 


( 2 . 21 ) 


Applying Eq. (2.18) to Eq. (2.12) gives the transformation functions for the velocity 
held: 


I'fi) 


{ 1. 

0 

w ) 


dZa CFz—H 

H 


\ dx H — Za ’ 

H-Za 



y w 


( 2 . 22 ) 


Using Eq. (2.20) and some algebraic manipulation the inverse transformation func¬ 
tions can be constructed: 


MH- 

x = x, y = y, z= [ -— 


fu] 


/ 

1, 

0 

w J 



dZa CTz—H 

H-Za 


v 

dx H ’ 

H 



/ 

V 



(2.23) 


(2.24) 
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C. SPATIAL DISCRETIZATION: CONTINUOUS GALERKIN 

Continuous Galerkin (CG) methods are the general family of methods that 
include: hnite element (FE), spectral element (SE), and spectral methods. CG meth¬ 
ods take complex geometries (elements) from physical space to computational space, 
using continuous basis functions that are used to approximate the solution to a given 
partial differential equation (PDE). To construct the problem relevant to this project, 
the generalized 2-D hyperbolic-elliptic PDE was first considered [7]: 

dq 2 

— + u-vq = z/V q 


where q = q{x,t), u = u{x), x = {x,zY, and v is the viscosity coefficient. Using 
Galerkin machinery, q and u were approximated using basis function expansions: 


Mm 

qN{x,t) = Y,'>Pji.x)qj{t) (2.25) 

Mm 

Fjv(T,t) = ^^7i(F)Mj(t) (2.26) 


where represents the number of points inside the quadrilateral elements (Mat = 
(A^ -|- 1)^), N is the order of the polynomial approximation, and ipj is the Lagrange 
polynomial basis functions. The approximations for q^ and un were then substituted 
into the PDE, multiplied by a test function, ipi, and integrated across the global 
domain, D, to get (weak integral form): 


/ + / ^j(m ■ Vgjv)dD = u ! 'ifjiW'^qNd^ VT G (2.27) 

Jn at Jn Jn 
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Instead of solving the global problem directly, 2-D local basis functions were con¬ 
structed. The 2-D local basis functions are dehned as: 


® hkiv) 


j,k = 0,1,...,N z = l,2,...,{N + lf 


where h is a 1-D local basis function and 0 is the tensor/outer product of the 1- 
D local basis functions. The 1-D Lagrange polynomial local basis function, using 
Legendre-Gauss-Lobatto interpolation points, was dehned by: 

n (1^) 

I = 0 

Additionally, in order to construct the basis functions and transition between physi¬ 
cal space {x, z) and computational space (.^, rj) requires knowledge of the metric terms: 

1 dz —1 dx 

dx \J\drj ’ dz \J\dr] 

dr] —1 dz dr] 1 dx 
dx \J\ d^ ’ dz \J\ d^ 

dx dz dx dz 
d^ dr] dr] d^ 

With Eq. (2.27) in local form it can now be converted using integration by parts 
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(IBP) giving: 


[ iJiV‘^qNdQ:e= [ V ■ (%l)iVqN)dVte - j V%l)i ■V{qN)dVte (2.28) 


then using the divergence theorem on Eq. (2.28) yields: 


f ipiV^qNdQe= f n ■ {ipiVqN)dre - f Vipi ■ V{qN)dQe (2.29) 

jQe JTe 


The result from Eq. (2.29) is then substituted back into the original PDE seen in Eq. 
(2.27), which produces: 



dqN 

dt 


dfle T 


%lJi{u- VqN)dQe 



{pJiVqN)dTe -ly ■V{qN)dile 

J 


Substituting the summation approximation for qjq and (Eqs. 2.25 and 2.26) yields: 



The resulting matrix problem is: 



where is the mass matrix, is the differentiation matrix (a discrete repre¬ 
sentation of V), is the boundary matrix, and is the Laplacian matrix. The 
discretized matrices are: 




(e) _ 


Mn 




1=1 


uiD. 


(e) 

ij 


Mm 

'^UJi 

1=1 


Mn 

'^KiUk 

k=l 




Q+l 

1=1 



Mn 


1=1 


where Q represents the quadrature used for evaluation of the integrals, which when 
using inexact integration is equal to N. Direct stiffness summation (DSS) was then 
used to construct the global problem: 


Ne 

Mjj = Am, 


(e) 


e=l 


LlJ — 


AL 


(e) 


e=l 


Dij = 


AD 


(e) 


e=l 


where Mu reduces to a diagonal matrix M/ when using inexact integration. 
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1. Equation Set 1 

Applying these operators to the Navier-Stokes Equation Set 1 (2.2-2.3) yields: 


^ = -ujMr^DijTTj - ^TTjMf^DjjUj 

—* 

= -iFjMj^DijUj - CpOMY^Djjnj - gMj^k 

^ = -uJMY^DijOj 


2. Equation Set 2 

Similarly to Set 1, the operators applied to the Navier-Stokes Equation Set 2 
(2.4-2.6) yield: 


dpi 

dt 


-Mj ^Dij{pu)j 


dui 


-uJMj ^DijUj - —Mj ^DijPj - gMj 
Pi 


ddi 


-ujMY^Dije,! 
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D. TEMPORAL DISCRETIZATION RUNGE-KUTTA 35 

Integrating nonlinear PDEs forward in time can lead to numerous problems 
(i.e. spurious oscillations, overshoots, progressive smearing, etc...) if the proper time 
integration scheme is not used [8]. For the proposed test cases, a strong stability pre¬ 
serving (SSP) Runge-Kutta (RK) method was selected and implemented as outlined 
in Ruuth and Spiteri [8]. SSP time integration methods have strong nonlinear stabil¬ 
ity properties, which make them optimal, for this thesis, for temporal discretization 
because of the nonlinearities present in the Euler and Navier-Stokes equations. The 
particular method used for time integration was an explicit third-order five-stage RK 
method (RK35). This method was chosen for its large stability region relative to 
other explicit methods of equal order. The RK35 method can be represented by: 


Q„ = Q‘"' 

7-1 

Qi = '^{oii^kQk) + / = l,2,...,s 

k=0 

gn+1 ^ q ( s ) 


where I is the stage (5 for RK35) and the coefficients aij and [3i are listed in the 
Appendix. 
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III. APPLIED COORDINATE 
TRANSFORMS 

A. EQUATION SET 1 

1. Perturbation Method 

This section will outline the expansion of the terms (application of the pertur¬ 
bation method) of set 1 of the Navier-Stokes Equations, where both tt and 6 will be 
split/decomposed into two components, the mean values (ft and 9) and their associ¬ 
ated perturbations (tt' and 9') such that: 


71 = 7l(z) + 7l'(x, Z, t) 


and 


9 = 9(z) -|- 9'(x, z, t) 


where the mean values satisfy a hydrostatically balanced atmosphere. After lineariza¬ 
tion, the pressure tendency Eq. (2.2) becomes: 


8(77 + 77 ) _ ,, R, ^ 

- - - - -I- M ■ V(7r -I- TT ) -I- (tt -I- TT )V ■ U = 0 

8t Cn, 


which, after simplihcation, becomes: 


877 ^ , 877 R, ^ 

—- h U ■ VtT -|- W— - 1 - (77 + TT )V ■ M = 0 

8t 8z Cv 


( 3 . 1 ) 
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Expansion of the terms of the momentum Eq. (2.2) yields: 


du 

'm 


+ u ■ Vu + Cp(0 + d')V{n + tt') 



which can be expanded to: 


du 

'm 




/ dn dn\ / dn' dn' \ 

y (9a; ’ (9^: J y (9a; ’ dz J 



from which, using the hydrostatic approximation ^ yields: 


du 

'm 


u • ^u Cp(^6 9'') 


g z- / dn' (97r'\ 



which can be simplihed to: 


du 

'm 


+ u ■ Vn + Cp{9 + 9') 


/ dn' (97r'\ 
ydx' dz J 


-gk-gjk 



and hnally becomes: 


du 

'm 


+ u ■ Vu + Cp{9 + 9')Vn' 



(3.2) 


Expansion of the terms of the thermodynamic energy Eq. (2.3) yields: 


d{9 + 9') 

m 


+ u-V{9 + 9') 


0 
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which further leads to: 


dQ' ^ de 

— + M ■ V9 + w— = 0 
at oz 


(3.3) 


2. Transform 

Using the basic machinery prescribed in Eqs. (2.13) - (2.22), the set of non¬ 
conservative Navier-Stokes Eq. (3.1 - 3.3) were transformed from x-z coordinates 
to x-az coordinates. The first machinery applied to the pressure tendency Eq. (3.1) 
was to change the vector dot products to transposes (i.e. u-V to fc^V), which yielded: 

— —h Vtt -|- w ——I- (tt + TT )V ■ u = 0 

at oz Cv 

substituting in Eq. (2.13) and Eq. (2.14) leads to: 

^ + x')(j^v)y j-‘s) = 0 

applying the linear algebra identity in Eq. (2.15) and the transformation for [w) 
yields: 


dn' 




( _dzs (Jz — H — Zs\ dn daz 

1^““^ H ^ ^ H ) daz dz 


+ ^{71+ 7r’){V^){u) = 0 
Cy 
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which simplihes to: 


dn' 


{uf {y)Ti' + 


( ^dzs(T^- H 
“9a; H 


+ w 


H-Zs \ 

H ) 


dn / H \ 

d(7z \H — Zg) 


+ —(tt + 7 r')(V)'^(-u) 


0 


and then to: 


dn' z, ~ , 

—+ «.Vx 


U 




H 


H-Zg 


dzg dn dn R, ~ 

— --h w- -h —(tt + tt )V ■ m = 0 

ox oaz oaz c,. 


and then further simplihed to: 


dn z, ~ , 

——h u ■ Vtt — u 
dt 


(Tz-H 

H-Zg 


dzg dn dn R, ~ 

— — + w- -^- (tt + tt )V ■ m = 0 

oaz dx Oaz c„ 


to hnally yield: 


dn z, ~ ^ dn R 

—— Vu-Vn + w— -1-(7r + 7 r)V-M = 0 (3.4) 

dt daz Cy 


Applying the same machinery as above to the momentum Eq. (3.2), the dot products 
were replaced, which yields: 


du 

'm 


+ ifVu + Cp{e + e')Vn' 



then u and V were replaced by their transforms (Eq. (2.13) and Eq. (2.14)) leading to: 

_ O! 

^ + {J-^'uf{J^V)u + Cp{e + e'){J^V)n' = gjk 
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applying the linear algebra identity in Eq. (2.15) yields: 


^ + (5)’’(J-')^(J^)(V)'S+ Cp(» + e’)(J'^Vy = gtk 


which simplihes to: 


at 

I\l _i 


— + ( v)w + c ,(0 + e'){j^ vy = g^k 


to hnally yield: 


— + u- Vu + Cp {0 + 6'')( j'^V)7r' = gjk 


Applying the machinery to thermodynamic energy, Eq. (3.3) yields: 


d9' de ^ 

— + yV9 + w— = 0 
ot oz 


substituting in Eq. (2.13) and Eq. (2.14) leads to: 


+ {j-yf {yV)9' + w 


09 Ooz 
Oz da. 


applying the linear algebra identity in Eq. (2.15) and the transformation for {w) 
yields: 


+ {uf {j-y {y){y)9'+ [-U 


.Oz. a. — H — z.\ 09 da 


Ox H 


H / Oa. Oz 
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which simplihes to: 


dt 


+ (€)^(V)0' + 


/ _dzs(T^- H 
( “9a; H 


+ w 


H-zA 

H ) 


89 [ H \ 
da^ \H — Zg) 


0 


and simplihes further, yielding: 


86' ~ ~ ^(Jz — H 8zs 89 

— + u.V9'-u^ -^ — 

8t H — Zg 8x 8az 


+ w 


89 

8az 


0 


and then to: 


89' ~ ~ , _az — H 8zg 89 ^ 89 

_ + -+u.— 


0 


to hnally yield: 


89' ~ ~ ^ 89 ^ 

— + u- V9 + w -— = 0 
8t 8az 


(3.6) 


3. Decomposition 

In order to discretize the governing equations and code them into Fortran 
90/95, the vector helds had to be decomposed into scalar components. The decom¬ 
position of the pressure tendency Eq. (3.4) became: 


97r' 

■ 97r' 97r'' 

8% R, 

8u 

8w 1 


“ o + 

8x 8az 

tc ^ -h in + n) 

8az Cy 

8x 

-1 

b 


= 0 
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(3.7) 



The momentum Eq. (3.5) decomposes into two separate equations (u and w), with 
the hrst equation [u) taking the form: 


du 

'm 


_du ^ du 
ox oaz 


+ Cp(0 + 9') 


dx 


/ dzs (Jz - H \ 

y (9a; H — Zs) 


dn' 

daz 


0 


(3.8) 


and the representation for [w) taking the form: 



dw 



+ w 


dw 

daz 


+ Cp(d + 9'^ 


\( ^ ] 


\H-zJ 

-1 

b 



(3.9) 


The decomposed thermodynamic energy Eq. (3.6) yields: 


'm 


d9' d9' 

dx daz 


+ w 


d9 

daz 


= 0 


(3.10) 


4. Application of the Galerkin Statement 

Using the Galerkin machinery outlined in the previous chapter, Eq. (3.7) - 
(3.10) were discretized yielding the continuous Galerkin Set 1 using yi-az coordinates 
(GGl x-az)-. 


Mn Qz^i, 

1=1 

Mn 
1=1 


Mn 

E 

k=l 


I (~ / I ~ / ~ - 


R 


(Mn 


-+ 


\fc=l 


d'ipj,i 

dx 


Uj + 


dipj^i 

daz 



(3.11) 
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1=1 


dt 


Mm 

\Ji\'4)i^i 

1=1 


^ I (~ , ~ 




dZs\ ^z^j H 
^ dx ) . H - Zsj 


da, ^ 


Qyj, 

1=1 


dt 


Mm 

'^Ul\Ji\ ili^i 


1=1 


/ f~ dipj^i . dipj^i 


da. 


^Mm 

-Cp I J2^k,l{^k + O'k) 


^k=l 


H 


^ ,j 


9^3,1 A , 
-^^j]+97f 


da. 


9, 


Mm QQ! 

Y.^1 = 

1=1 ^ 


Mm 

'^UJl \Jl\iJi,l 


1=1 


Mm 


^Ha , -r. ^Ha'\ . 7 . S’l’iJ 


EV'« + 


k=l 


dx ^ da. ■'! da. '’ 


The existing CGI x-z code (in Fortran 90/95) [1] was modified rising Eqs. 
(3.14) and then used for the test cases outlined in the next chapter. 


(3.12) 

(3.13) 

(3.14) 

(3.11) - 



B. EQUATION SET 2 


1. Perturbation Method 

This section will outline the expansion of the terms of set 2 of the Navier- 
Stokes Equations, where (as seen in set 1) both p and 9 will be split/decomposed into 
two components, the reference values (p and 6) and their associated perturbations (p' 
and 9'). As seen in set 1: 


p = p{z) + p\x,z,t) 


and 


9 = 9(z) + 9'(x, z, t) 


where the reference values again satisfy a hydrostatically balanced atmosphere. In 
order to linearize, the mass Eq. (2.4) the product rule is first applied, yielding: 


dp ^ ^ vv- 

M ■ Vp + pV ■ M 
at 


0 


followed by expansion of the terms, becomes: 

^k^ + ^.V(p + p') + (p + p')V-M = 0 


which, after simplihcation, yields: 
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(3.15) 


% + u-Vp' + w^ + {p + p')V ■ S = 0 

at oz 


By setting fi = 0, the momentum Eq. (2.6) becomes: 


du r 

— + u- Vu + -VP = —gk 
at p 


and followed by an expansion of the terms, yields: 


1 — 

— + U-VU+- --V(P + P0 

dt {p + p') 


which leads to: 


du^ 1 , 1 dP - 


from which, using the hydrostatic approximation ^ = —pg yields: 


du ^ 1 „ PQ r r 


which, after simplihcation, becomes: 


^ M ■ Vm + ^ VP' + k = 0 

dt {p + p') (p + p') 
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and finally: 


du ^ ^ 

— + M ■ Vm + 

ot 


1 

(p + p') 


VP' 


p'g 

{p + p') 


k 


By setting fi = 0, the thermodynamic energy Eq. (2.6) becomes: 


^ + VP = 0 
ot 


and followed by an expansion of the terms, yields: 


d{e + 6') 


+ u-v{e + e') 


0 


(3.16) 


which further leads to: 


de' ^ 89 ^ 

— + u- V6 + w— = 0 
at dz 


(3.17) 


2. Transform 

Using the basic machinery described in Eq. (2.13) - (2.22), the set of non¬ 
conservative Navier-Stokes Eq. (3.15 - 3.17) were transformed from x-z coordinates 
to x-(Tz coordinates. The hrst machinery applied to the pressure tendency Eq. (3.15) 
was to change the vector dot products to transposes (i.e. u-V to n^V), which yielded: 

^ -h iFVp' + + (p + p’)V'^u = 0 

at dz 
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substituting in Eq. (2.13) and Eq. (2.14) leads to: 


dp' 




dp daz 
dz daz 


+ {p + p'){J^VfJ-^^ 


0 


applying the linear algebra identity in Eq. (2.15) and the transformation for [w) 
yields: 


dp' 




( ^dz^Oz-H 

“Sa; H 


+ w 


H-Zs \ 

H J 


dp daz 
daz dz 


+ {p + p'){yf J-^~u = Q 


which simplihes to: 


dp' 


+ {^fVp + 


^ dzs az — H — Zs\ dp ( H \ 

H ^ ^ H ) ^ \H-Zs) 


+ (p + p')(Vf'u 


0 


to hnally yield: 


^ + Vp' + «;^ + (p + p')V-^=0 
dt da. 


(3.18) 


Applying the same machinery as above to the momentum Eq. (3.16), the dot prod¬ 
ucts were replaced, which yields: 


du 

dt 


w Vu + 


(p + p'\ 


-VP' 


p'g 

(p + p') 


k 
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then u and V were replaced by their transforms (Eq. (2.13) and Eq. (2.14)) leading to: 


at ip + p') (p + p') 


applying the linear algebra identity in Eq. (2.15) yields: 


^ t- 1\T jT^rt , 1 P'9 


^+u{j-yrvu + -^—rvp'= , 

dt (p + pO (p + pO 


which simplihes to: 


^ + / w+ k 

dt (p + pO (p + p') 


to hnally yield: 


du z, 1 p'q r 

— + u - Vu + — - -J^VP = - - -k 

dt (p + pO (p + p') 


Applying the machinery to thermodynamic energy Eq. (3.17) yields: 


d9' de 

— + u^V6 + w— = 0 
dt dz 


substituting in Eq. (2.13) and Eq. (2.14) leads to: 


dO' 

'm 


. _ 1 Tlx 7^ /y ~ . QO dCT^ 

+ (J ^u)^J^V9 + w—t :— = 0 

dz da. 


(3.19) 


33 



applying the linear algebra identity in Eq. (2.15) and the transformation for [w) 
yields: 


86' 


(ny j^ve' + 


/ ^dzsCTz- H 

( “ax H 


+ w 


H-zA 

H ) 


89 8(7z 
8az 8z 


0 


which simplihes to: 


89' 


{JiyV9' 


^8zsaz-H .H-Zs\ 89 ( H \ 
H ^ ^ H ) ^ \H-Zs) 


0 


to hnally yield: 


90' r. ~ 
+ u ■ V9 + w 
8t 


89 

8az 


0 


(3.20) 


3. Decomposition 

Similar to CGI, in order to discretize the governing equations and code them 
into Fortran 90/95, they had to be decomposed into scalar components. The decom¬ 
position of the mass Eq. (3.18) became: 


8p' \^8p' ^ 8p' 


-|- w 


8p 

8az 


+ (P + P') 


8u 8w 
8x ^ 8a z 


0 


(3.21) 


The momentum Eq. (3.19) decomposes into two separate equations {u and w), with 
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the first equation [u) taking the form: 


du 


^du ^ du 
ox oaz 


+ 


(p + P') 


dP' ( dzsaz-H \ dP' 
dx ^ \ dx P[ — z.) da. 


= 0 (3.22) 


and the representation for (w) taking the form: 


p'a 

{p + P') 

The decomposed thermodynamic energy Eq. (3.20) yields: 


dw 


^ dw _ dw 
dx da. 


+ 


(p + P') 


H 


H — z.J da 


dP' 


(3.23) 


dt 


■ dO' 
u ——h w 
dx 


de' 

daz 


+ w 


de 

da. 


0 


(3.24) 


4. Application of the Galerkin Statement 

Using the Galerkin machinery outlined in the previous chapter, Eqs. (3.21) - 
(3.24) were discretized yielding the continuous Galerkin Set 2 using yi-az coordinates 
(GG2 x-a^): 


Mn QJ, 

1=1 ^ 


Mm 

1=1 


Mn 

I “fc 


k=l 


d'dj^i 

dx 


+ 


' dz. 


a 




dx J . H 


H\ dip,,I 

da. 




Pi 


( H \ dipj,i , ^ dipj,i _ " 


/Mm 


Y.^kAPk +P'k)\ 


<k=l 


dipj I ^ dipj I _ ' 
-i-hin . j_. 

dx ^ da. ^, 


(3.25) 
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The existing CG2 x-z code (in Fortran 90/95) was modified using Eq. (3.25) - (3.28) 
and then used for the test cases outlined in the next chapter. 
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IV. 


TEST CASES 


A. OVERVIEW 

Three specific cases where chosen to evaluate the a-coordinate transforma¬ 
tions: the rising thermal bubble, the linear hydrostatic mountain, and the linear non¬ 
hydrostatic mountain. The rising thermal bubble, case 1, was chosen as a benchmark 
in order to insure that in the absence of terrain the model dynamics still functioned 
properly (i.e. the a-coordinate transformed equations reduce to the original govern¬ 
ing equations) yielding the same results as the unmodified source code [1]. The linear 
hydrostatic mountain, case 2, and the linear non-hydrostatic mountain, case 3, were 
chosen to evaluate the advantages and disadvantages of the x-a governing equation 
in relation to their x-z formulation in both hydrostatic and non-hydrostatic environ¬ 
ments using the same terrain. This chapter will outline the test case assumptions, 
additional information derived for the cr-coordinates, and previous results that will 
be used for comparison. 

B. CASE 1: RISING THERMAL BUBBLE 

This test case represents a highly nonlinear non-hydrostatic flow problem, by 
demonstrating the evolution of an initially at rest warm air bubble, relative to the 
surrounding environment, in a hydrostatic constant potential temperature environ¬ 
ment as resolved by the Navier-Stokes equations. As the warm air bubble rises, it will 
deform, taking the shape of a mushroom cloud, due to the shearing motion caused by 
the resulting differential vertical velocity field. The initial distribution of potential 
temperature perturbations in order to drive vertical motion is defined by: 


e' 


0 

Be 

2 


1 -|- COS 



r > Tc, 

r < Tc, 
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where 9c = 0.5°C, tTc is the trigonometric constant, r = \J{x — XcY + {z — Zc^ with 
the constants: 9 = 300 K, Vc = 250 m, and {xc, Zc) = (500,350) m. [1] The domain 
of interest is (x,z) G [0,1000]^ m integrated over the time interval, t G [0, 700] s. 
Additionally, the boundary conditions for all boundaries are no-flux. 

The utility of running this case is to insure that the a transformed governing 
equations yield the same results as the x-z set of the governing equations when the 
terrain is flat: 


Zsurf 0 


where Zsurf is the height of the surface. The surface heights are vital in deriving the 
(T-coordinates (see equation 2.16). The slope of the surface must also be calculated 
for the transformation, which in this case leads to: 

9zg^cf _p, 

dx 

With a flat surface in combination with a slope of 0, the transformation of the coor¬ 
dinate system (x-Za) reduces to x-z, making both the modihed and unmodihed codes 
yield identical results. 

The unmodihed source code was run in order to establish a baseline to evaluate 
the transformed source code. Three resolutions were used: 20 m (2601 grid points), 10 
m (10201 grid points), and 5 m (40401 grid points). The resulting data for potential 
temperature perturbations was then plotted for each resolution after 700 s of model 
integration (see Figures 2 and 3). As seen in Figures 2 and 3, the warm air bubble 
did deform into a mushroom type formation, while maintaining symmetry. 
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Figure 2. Case 1: Rising Thermal Bubble. Resulting potential temperature pertur¬ 
bations using unmodified CGI source codes [1] after 700 s for resolutions: (a) 20, (b) 
10, and (c) 5 m. All cases were run using lOtli order polynomials, with contours from 
-0.05 to 0.525 K with an interval of 0.025 K. 
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Figure 3. Case 1: Rising Thermal Bubble. Resulting potential temperature pertur¬ 
bations using unmodified CG2 source codes [1] after 700 s for resolutions: (a) 20, (b) 
10, and (c) 5 m. All cases were run using lOtli order polynomials, with contours from 
-0.05 to 0.525 K with an interval of 0.025 K. 
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C. CASE 2: LINEAR HYDROSTATIC MOUNTAIN 

This test case was chosen in order to evaluate the model performance when 
a simple terrain feature is introduced in a linear hydrostatic flow environment with 
steady inflow and outflow lateral boundary conditions. Over time, a steady-state 
mountain wave over a single peak should form if the model dynamics are resolving 
the scenario accurately. Initially, the constant mean horizontal velocity is set io u = 
20 m/s, the mean vertical velocity is set to tc = 0 m/s, and the atmosphere is set 
to isothermal with a constant mean temperature of T = 250 K.[l] Since this case 
is isothermal, the buoyancy frequency or Brunt-Vdisala frequency A/” ^ = g-^{ln9) 
reduces to A/" = which simplihes the Exner pressure to: 


Q 

7f = e 


The domain of interest is (x,z) G [0, 240,000] x [0, 30,000] m integrated over the time 
interval, t G [0,10] h. For this test case, the terrain (the versiera di Agnesi mountain 
prohle) is represented by: 


he 

(i + (^)) 

where Zgurf is the height of the surface, and Xc, and Oc are constants {he = 1 m, Xc 
= 120,000 m, and Oc = 10,000 m). Additionally, the bottom boundary conditions are 
no-flux, while the top and lateral boundaries use a non-reflecting boundary condition 
[1]. For further information on such boundary conditions, see the papers by Durran 
and Klemp [9], Giraldo and Restelli [1], and Dea, Giraldo, and Neta [10]. However 
non-reflecting boundary conditions are outside the scope of this thesis and will not be 
mentioned further. Unlike the previous case, the surface heights do vary with x and 
are even more vital in deriving the a-coordinates. The surface plot can be seen in 
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Figure 4.a, on which it can be seen (noticing the z-axis scale) that the peak is fairly 
low and has smooth geometry. Since the surface height does vary with x, the slope of 
the surface must also be calculated for the transformation, which in this case leads to: 


^surf he j 1 ~\~ 


X — XC 


ac 


-1 


^surf he 


(ae)^ + {x^ — 2{xc){x) + (xe)^) 


ae 


-1 


Zsurf = {hc){acy (acY + x^ — 2{xc){x) + (xc) 


n-i 


dzsurf _ {-l){hc){acY{2x - 2{xc)) 

dx [{acY + x‘^ — 2{xc){x) + {xcYY 

dzsurf {—2){hc){acY{x — xc) 
dx [(ac)2 + {x — xc)‘^Y 

where the surface slope (—|^) plot can be seen in Figure 4.b, on which it can be 
seen (noticing the z-axis scale) that the slope of the peak is fairly flat and again has 
smooth geometry. 

The unmodihed source code was run (colored lines) and compared with the 
analytic solutions (black lines), generated for the unmodihed source code and de¬ 
veloped by Giraldo and Restelli [1], in order to establish a baseline to evaluate the 
transformed source code. One resolution was used: 1200 m (in x) and 240 m (in z) 
(20301 grid points). The resulting data for vertical and horizontal velocity was then 
plotted for the resolution after 10 h (see Figure 5). As seen in Figure 5, a steady-state 
mountain wave did form in the appropriate region and will be compared against the 
transformed governing equation sets. 
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Figure 5. Case 2: Linear Hydrostatic Mountain. Resulting horizontal velocity (a) and 
vertical velocity (b) using unmodified (i) CGI and (ii) CG2 source codes [1] after 10 
h for the resolution of 1200 m (in x) and 240 m (in z) (colored lines) plotted with the 
analytic solution (black lines). All cases were run using 10th order polynomials, with 
contours from -0.025 to 0.025 ms“^ with an interval of 0.005 ms“^, (a), and -0.005 to 
0.005 ms“^ with an interval of 0.0005 ms“^, (b). 
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D. CASE 3: LINEAR NON-HYDRO STATIC MOUNTAIN 

The final test case was chosen in order to evaluate the model performance when 
a simple terrain feature is introduced in a linear non-hydrostatic flow environment 
with steady inflow and outflow lateral boundary conditions. Over time, a steady-state 
mountain wave over a single peak should form if the model dynamics are resolving the 
scenario accurately. Initially, the constant mean horizontal velocity is set to h = 10 
m/s, the mean vertical velocity is set to ta = 0 m/s, and the atmosphere is uniformly 
stratified with a Brunt-Vaisalci frequency of Af = 0.01 /s. This test case uses the 
same terrain as case 2 (see Figure 4). The constants used for this case were he = f 
m, Xc = 72,000 m, Oc = 1,000 m, and 9o = 280K. The domain of interest is (x,z) G [0, 
144,000] X [0, 30,000] m integrated over the time interval, t G [0, 5] h. Like case 2, 
the bottom boundary conditions are no-flux and the top and lateral boundaries are 
non-reflecting. [1] 

The unmodified source code was again run (colored lines) and compared with 
the analytic solution [1] (black lines) in order to establish a baseline to evaluate the 
transformed source code. One resolution was used: 360 m (in x) and 300 m (in z) 
(40501 grid points). The resulting data for vertical and horizontal velocity was then 
plotted after 5 h (see Figure 6). As seen in Figure 6, a steady-state mountain wave 
did form in the appropriate region and will be compared against the transformed 
governing equation sets. 
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Figure 6. Case 3: Linear Non-Hydrostatic Mountain. Resulting horizontal velocity 
(a) and vertical velocity (b) using unmodified (i) CGI and (ii) CG2 source codes [1] 
after 5 h for the resolution of 360 m (in x) and 300 m (in z) (colored lines) plotted with 
the analytic solution (black lines). All cases were run using 10th order polynomials, 
with contours from -0.025 to 0.025 ms“^ with an interval of 0.005 ms“^, (a), and 
-0.005 to 0.005 ms“^ with an interval of 0.0005 ms“^, (b). 
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V. 


RESULTS 


A. CASE 1: RISING THERMAL BUBBLE 

1. Accuracy and Comparison 

The modified source code (containing the a-coordinate transformed governing 
equations) was run using the same three resolutions that were used for the unmodified 
code: 20 m (2601 grid points), 10 m (10201 grid points), and 5 m (40401 grid points). 
The resulting data for potential temperature was then plotted for each resolution after 
integrating forward 700 s (see Figure 7). As seen in Figure 7, the warm air bubble 
did deform into a mushroom type formation, while maintaining symmetry, appearing 
to duplicate the results seen by the x-z code (see figure 2). 

To further evaluate the similarities in the x-z and x-a^ model solutions nine 
variables (using 5 m resolution) were compared: the maximum Exner pressure pertur¬ 
bation (7r(„£ja, (unitless)), the minimum Exner pressure perturbation (ttA^ (unitless)), 
the maximum horizontal wind velocity {umax (ins”^)), the minimum horizontal wind 
velocity {umin (ms“^)), the maximum vertical wind velocity {wmax (ms“^)), the mini¬ 
mum vertical wind velocity {wmin (^is”^)), the maximum potential temperature per¬ 
turbation (K)), the minimum potential temperature perturbation (6*A„ (K)), 

and CPU time (s). These values were compare to better observe the relative maxi¬ 
mum and minimum values occurring for each model run. In order to have consistent 
CPU constraints both x-z and x-a^ for set 1 code were run in parallel on the same 
machine (NPS Math Department 32 Processor Apple Cluster: Riemann) and started 
at the same time, and then repeated for set 2. The resulting values can be seen in Ta¬ 
ble I. All four sets of code converge to nearly identical solutions. The only difference 
are the following: minor differences in the 7r[nax 1 2; and differences in 

0'^^^ and between x-z coordinates and x-a^ coordinates for both sets. The major 
difference between the models is the CPU time (in seconds), which was consistently 
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slower for the x-a^ ( 28.24 % increase in time for set 1 and 15.58 % increase in time 
for set 2) relative to their x-z counterparts. This difference is caused by the additional 
x-az data being passed between program functions and additional calculations being 
performed ( 0(5 * [Nt] * [Ne] * * [8 * (iV + 1) + 26])), where Nt is the number 

of time steps, Ne is the number of elements, and Mjv is the number of quadrature 
points or order iV + 1). For the 5 m resolution case, there were 5.2101x10^® additional 
operations. The horizontal velocities Umax and Umin also indicate that symmetry is 
maintained for all four model sets. 

To further evaluate the performance of the four model runs, the Root Mean 
Squared Error (RMSE) was determined for each run in relation to the analytic solu¬ 
tion (see Table II) for the four variables: tt, u, w, and 9. The RMSE was calculated 
using: 


\RMS 


\ 


1 Np 

_ ^ ^ numerical _ q analytic'^^ 

i=i 


where q represents the given state variable vector, Np is the number of points, and the 
analytic solution for this case was dehned by the unmodihed source code numerical 


Table I. Case 1: Rising Thermal Bubble. Comparison of modihed and unmodihed 


CGI and CG2 using 5 m resolution after 700s. 


Model 

CGI x-z 

CGI x-az 

CG2 x-z 

CG2 x-a. 

TT^ 

''max 

0.9364x10"^ 

0.9364x10“^ 

0.9365x10-® 

0.9365x10-® 

Tl' ■ 
mm 

-0.1195x10-^ 

-0.1195x10-^ 

-0.1195x10-^ 

-0.1195x10-^ 

Umax (mS“^) 

0.2079x10^ 

0.2079x101 

0.2079x101 

0.2079x101 

Umin (mS“^) 

-0.2079x10^ 

-0.2079x101 

-0.2079x101 

-0.2079x101 

Wmax (mS“^) 

0.2536x10^ 

0.2536x101 

0.2536x101 

0.2536x101 

Wmin (mS"^) 

-0.1912x10^ 

-0.1912x101 

-0.1912x101 

-0.1912x101 

O'max (K) 

0.5713x10° 

0.5715x10° 

0.5713x10° 

0.5715x10° 

OLn (K) 

-0.9736x10-1 

-0.9735x10-1 

-0.9736x10-1 

-0.9735x10-1 

CPU time (s) 

0.1027x10^ 

0.1317x10® 

0.1322x10® 

0.1528x10® 
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Table II. Case 1: Rising Thermal Bubble RMSE. Root-mean-squared errors for the 
four primary state variables, for the modihed codes using the unmodihed code so¬ 
lutions as the analytic solutions, after 700 s using 5 m resolution and 10th order 
polynomials. 


Model 

n 

u (ms 1) 

w (ms 1) 

0(K) 

CGI x-az 
CG2 x-az 

2.5931x10-1^ 

2.2034x10-1® 

9.2434x10-1® 

5.2908x10-12 

1.6068x10-® 

8.5529x10-12 

1.8639x10-® 

9.6233x10-12 


solutions in order to better compare x-z solution to the x-a^ solution. Both sets 1 
and 2 exhibited a smaller error for all four variables of interest, with set 2 having the 
lower RMSE values overall. 

In addition to comparing extremes of the selected helds, a one-dimensional 
vertical temperature profile along the centerline (x = 500 m) was plotted for each 
resolution to better discern the differences between the various models (see Figure 
9). As the resolution increased from 20 m to 5 m, there is a noticeable increase in 
the sharpness of the temperature distributions. As expected, for all three resolutions, 
there is no significant difference between the temperatures along the centerline for 
each model at that specihc resolution. 

2. Conclusions 

All four model runs for the rising thermal bubble did deform into a mushroom 
type formation, while maintaining symmetry (as indicated by the distribution of hor¬ 
izontal velocities). Additionally, the models converged to nearly identical solutions, 
as demonstrated by the min and max values and the RMSE. Overall, for both sets 1 
and 2 of the Navier-Stokes equations, the x-az coordinates do reduce to x-z, showing 
that in the absence of terrain that the model dynamic still function properly, but did 
show a noticeable increase in computational expense. These results validate the first 
stage of the code evaluation. 
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Figure 7. Case 1: Rising Thermal Bubble. Resulting potential temperature pertur¬ 
bations using modified CGI source codes [1] after 700 s for resolutions: (a) 20, (b) 
10, and (c) 5 m. All cases were run using 10th order polynomials, with contours from 
-0.05 to 0.525 K with an interval of 0.025 K. 
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Figure 8. Case 1: Rising Thermal Bubble. Resulting potential temperature pertur¬ 
bations using modified CG2 source codes [1] after 700 s for resolutions: (a) 20, (b) 
10, and (c) 5 m. All cases were run using 10th order polynomials, with contours from 
-0.05 to 0.525 K with an interval of 0.025 K. 
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Figure 9. Case 1: Rising Thermal Bubble. Resulting potential temperature perturba¬ 
tions for all four models along the vertical axis (x = 500 m) after 700 s for resolutions: 
(a) 20, (b) 10, and (c) 5 m. All cases were run using 10th order polynomials. 
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B. CASE 2: LINEAR HYDROSTATIC MOUNTAIN 


1. Accuracy and Comparison 

The modified source code was run and the numerical solution after 10 hours 
(colored lines) was plotted with and compared to the analytic solution [1] (black lines) 
(see Figure 10). Only one resolution was used: 1200 m (in x) and 240 m (in z) (20301 
grid points). The resulting data for the vertical and horizontal velocity was then used 
for the comparison (see Table III). As seen in Figure 10, a steady-state mountain 
wave did form in the appropriate region closely modeling the analytic solution. 

Visually comparing Figure 5 to Figure 10, indicated fairly identical solutions, 
and to further evaluate the similarities in the x-z and x-cTz model solutions nine 
variables were again compared: ttV^, u^ax, Umin, Wmax, Wmin, d'max^ and 

CPU time (in seconds). Similar to case 1, both x-z and x-a^ for set 1 code were run 
in parallel and started at the same time, and then repeated for set 2. The resulting 
values can be seen in Table III. All four sets of code converge to similar solutions, with 
only minor differences overall, and the most predominate differences in the Wmax and 
Wmin between x-z coordinates and x-Uz coordinates, where x-z had greater magnitude 
values for Wmin and x-Uz had greater magnitude values for Wmax- As seen in case 1, 
another major difference between the models is the CPU time, which was consistently 
slower for the x-Uz ( 30.03 % increase in time for set 1 and 15.67 % increase in time 
for set 2) relative to their x-z counterparts. This difference is again being caused 
by the additional x-Uz data being passed between program functions and additional 
calculations being performed ( 0(5 * [Nt] * [V^] * * [8 * (V -|- 1) -|- 26])). For 

this case, there were 1.8793x10^^ additional operations. 

To further evaluate the performance of the four model runs, the RMSE was 
determined for each run in relation to the analytic solution (see Table IV) for the four 
variables: n, u, w, and 9. The analytic solution was interpolated to the model grid 
using a cubic approximation. Once the domain of interest was defined and the nu- 
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Table III. Case 2: Linear Hydrostatic Mountain. Comparison of modified and un¬ 
modified CGI and CG2 using 1200 m (in x) and 240 m (in z) resolution after 10 
hours._ 


Model 

CGI x-z 

CGI x-dx 

CG2 x-z 

CG2 x-az 

TT^ 

''max 

0.1792x10-^ 

0.1786x10-® 

0.1781x10-® 

0.1785x10-® 

TT^ ■ 

''min 

-0.1859x10-® 

-0.1859x10-® 

-0.1840x10-® 

-0.1846x10-® 

^max (rilS ) 

0.4072x10-^ 

0.4070x10-1 

0.4034x10-1 

0.4044x10-1 

Umin (mS"^) 

-0.3506x10-1 

-0.3484x10-1 

-0.3474x10-1 

-0.3473x10-1 

Wmax (mS“^) 

0.4229x10-2 

0.4656x10-2 

0.4245x10-2 

0.4686x10-2 

Wmin (mS"^) 

-0.5199x10-2 

-0.4279x10-2 

-0.5214x10-2 

-0.4313x10-2 

O'max (K) 

0.2400x10-1 

0.2386x10-1 

0.2380x10-1 

0.2382x10-1 

OLn (K) 

-0.3184x10-1 

-0.3124x10-1 

-0.3157x10-1 

-0.3115x10-1 

CPU time (s) 

0.4519x10^ 

0.5876x10^ 

0.6463x10^ 

0.7476x10^ 


merical and analytic solutions were on matching grids a bootstrap (random sampling) 
method was used to calculate the 95 % conhdence interval (Cl) [11]. The 95% Cl 
was needed in order to determine if the differences observed in the RMSE values were 
indeed significant. The bootstrap method built the 95% Cl by creating a domain (of 
equal size to the domain of interest) that was populated with random samples from 
the original numerical and analytic solution pairs, and then the RMSE for the new 
domain was calculated and stored. This process was iterated 10,000 times, storing 
the RMSE from each iteration. The derived RMSE values were then sorted and the 
values at 2.5% and 97.5% of the distribution were taken and compared to the original 
RMSE in order to establish the 95% CL The resulting conhdence intervals indicate 
that the differences between the state variable RMS errors are signihcant, since the 
ranges do not overlap. For set 1, x-z coordinates exhibited a signihcantly smaller 
error for all four variables of interest. For set 2, x-a^ had lower RMSE values for n 
and u, but larger RMSE values for w and 6. Taking a closer look at the w RMSE 
with respect to Wmax shows that CGI x-z has a 1.32% error as compared to 4.24% 
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Table IV. Case 2: Linear Hydrostatic Mountain RMSE. Root-mean-squared errors for 
the four primary state variables, for both the modihed and unmodihed codes, after 
10 h using 1200 m (in x) and 240 m (in z) resolution and 10th order polynomials. 


Model 

n 

95% Gonhdence Interval 

CGI x-z 

1.2263x10-^ 

+ 8.9354x10-1'^ 

- 8.9757x10-1° 

CGI x-az 

1.2898x10-^ 

+ 9.5565x10-1° 

- 9.9149x10-1° 

CG2 x-z 

1.2121x10-^ 

+ 9.1186x10-1° 

- 9.1327x10-1° 

CG2 x-az 

1.2200x10-^ 

+ 9.0255x10-1° 

- 9.1227x10-1° 

Model 

u (ms“^) 

95% Gonhdence Interval 

CGI x-z 

2.2561x10"^ 

+ 1.7930x10-3 

- 1.8157x10-3 

CGI x-az 

2.4089x10-3 

+ 1.9437x10-3 

- 1.9878x10-3 

CG2 x-z 

2.2481x10-3 

+ 1.8488x10-3 

- 1.8337x10-3 

CG2 x-az 

2.2808x10-3 

+ 1.8327x10-3 

- 1.8985x10-3 

Model 

w (ms-^) 

95% Gonhdence Interval 

CGI x-z 

6.8468x10-3 

+ 6.6744x10-^ 

- 6.6389x10-^ 

CGI x-az 

1.9757x10-^ 

+ 2.1074x10-3 

- 2.1468x10-3 

CG2 x-z 

5.1424x10-3 

+ 5.0225x10-^ 

- 5.0915x10-^ 

CG2 x-a. 

1.9311x10-^ 

+ 2.0131x10-3 

- 2.0667x10-3 

Model 

e{K) 

95% Gonhdence Interval 

GGl x-z 

1.6565x10-3 

+ 1.3026x10-3 

- 1.3132x10-3 

GGl x-az 

1.8506x10-3 

+1.4243x10-3 

- 1.4748x10-3 

GG2 x-z 

1.6360x10-3 

+ 1.4812x10-3 

- 1.4410x10-3 

GG2 x-az 

1.7208x10-3 

+ 1.4790x10-3 

- 1.5027x10-3 


for CGI x-(T^ and that CG2 x-z has a 0.99% error as compared to 4.12% for CG2 
x-a^, which is substantial for both sets. 

In order to verify that a steady-state solution for the linear hydrostatic moun¬ 
tain case was achieved, the momentum flux was derived using: 


OO 

m{z) = j p{z)u{x, z)w{x, z)dx (5.1) 

— OO 


where p{z) is the reference density as a function of height {z). [12] Additionally, the 
analytic hydrostatic momentum flux (from linear theory) is given by: 
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(5.2) 



4 


PsUsMhl 


where denotes the analytic hydrostatic momentum flux, ps is the reference density 
at the surface, Us is the horizontal velocity values at the surface, M is the Brunt- 
Vaiscila frequency, and he is the height of the mountain. [1] The momentum flux was 
then normalized by m{z)/m^{z), and will be the value discussed in this thesis. In 
Figure 11, the normalized momentum flux for all four model runs, using 1200 m (in 
x) and 240 m (in z) resolution, was plotted for 2 h, 4 h, 6 h, 8 h, and 10 h. As seen in 
Figure 11, all four model runs converge to a steady state solution after 10 h. Of the 
four model simulations, the CGI x-z and CG2 x-z models yielded results that were 
far better than the x-az models, with values between 0.95 and 1.01 versus between 
0.65 and 1.15 as seen for the x-a^ runs. 

2. Conclusions 

All four model runs for the linear hydrostatic mountain case did develop a 
steady-state mountain wave over a single peak indicating that the model dynamics 
are resolving the scenario accurately. Additionally, the models converged to nearly 
identical solutions. Overall, for both set 1 and set 2 of the Navier-Stokes equations, 
the x-az coordinates performed slightly worse than their x-z counterparts (reflected 
in the RMSE, the normalized momentum flux, and most pronounced in w with an 
approximately four times larger RMSE with respect to Wmax ) and as seen in the 
previous case did show a noticeable increase in computational expense. These results 
appear to indicate that when the model dynamics are resolved using purely explicit 
time integration methods that there is no clear advantage in using x-az over x-z and 
will be further explored into the next case. 
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Figure 10. Case 2: Linear Hydrostatic Mountain. Resulting horizontal velocity (a) 
and vertical velocity (b) using modified (i) CGI and (ii) CG2 source codes [1] after 10 
h for the resolution of 1200 m (in x) and 240 m (in z) (colored lines) plotted with the 
analytic solution (black lines). All cases were run using 10th order polynomials, with 
contours from -0.025 to 0.025 ms“^ with an interval of 0.005 ms“^, (a), and -0.005 to 
0.005 ms“^ with an interval of 0.0005 ms“^, (b). 
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i) 


ii) 




iii) iv) 

Figure 11. Case 2: Linear Hydrostatic Mountain. Normalized momentum flux for 
the resolution of 1200 m (in x) and 240 m (in z), at times 2 h, 4 h, 6 h, 8 h, and 10 h 
for the four model runs: (i) CGI x-z, (ii) CGI x-cxz, (iii) CG2 x-z, and (iv) CG2 x-a^. 
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C. CASE 3: LINEAR NON-HYDROSTATIC MOUNTAIN 


1. Accuracy and Comparison 

For the final case, both CGI x-a^ and CG2 x-a^ were run and the numerical 
solution after 5 hours (colored lines) was plotted and compared to the analytic solution 
[1] (black lines) (see Figure 12). Similar to case 2, only one resolution was used: 360 
m (in x) and 300 m (in z) (40501 grid points). The resulting data for the vertical 
and horizontal velocity was then used for the comparison (see Table V). As seen in 
Figure 12, a steady-state mountain wave did form in the appropriate region, there 
were some apparent contrasts with the analytic solution for CGI x-a^ and CG2 x-a^. 

Visually comparing Figure 12.i.a,b to Figure 12.ii.a,b indicated fairly similar 
solutions for the horizontal velocities (a). Comparing the vertical velocities (b) re¬ 
vealed a similar pattern in the placement of the steady state waves between x-z and 
x-a^, but the x-a^ coordinates appear to have significantly stronger oscillations than 
the x-z coordinates. To further evaluate the x-z and x-cr^ model solutions nine vari¬ 
ables were again compared: TtV^, Umax, Umin, Wmax, Wmin, O'max^ and CPU 

time (in seconds). Similar to case 1 and case 2, both x-z and x-a^ for set 1 code 
were run in parallel and started at the same time, and then repeated for set 2. The 
resulting values can be seen in Table V. All four sets of code converge to steady state 
solutions, with larger differences than what were observed in case 2, but had only 
minor differences in T^'max^ Umax, and Umin- The most predominate differences 

were observed in the Wmax and Wmin between x-z coordinates and x-a^ coordinates, 
where x-a^ had greater magnitude values for both state variables. There was also a 
noticeable difference in 9'max and between x-z coordinates and x-a^ coordinates 
for both sets. Consistent with the previous two cases, the CPU time was slower for 
both x-az model simulations ( 52.06 % increase in time for set 1 and 15.12 % increase 
in time for set 2) relative to the x-z simulations. This difference is again being caused 
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by the additional x-az data being passed between program fnnctions and additional 
calculations being performed ( 0(5 * [Nt] * [iV^;] * * [8 * (iV + 1) + 26])). For 

this case, there were 1.8700x10^® additional operations. 



Figure 12. Case 3: Linear Non-Hydrostatic Mountain. Resulting horizontal velocity 
(a) and vertical velocity (b) using both modified, (ii), and unmodified, (i), CGI source 
codes [1] after 5 h for the resolution of 360 m (in x) and 300 m (in z) (colored lines) 
plotted with the analytic solution (black lines). All cases were run using 10th order 
polynomials, with contours from -0.025 to 0.025 ms“^ with an interval of 0.005 ms“^, 
(a), and -0.005 to 0.005 ms“^ with an interval of 0.0005 ms“^, (b). 
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Figure 13. Case 3: Linear Non-Hydrostatic Mountain. Resulting horizontal velocity 
(a) and vertical velocity (b) using both modified, (ii), and unmodified, (i), CG2 source 
codes [1] after 5 h for the resolution of 360 m (in x) and 300 m (in z) (colored lines) 
plotted with the analytic solution (black lines). All cases were run using 10th order 
polynomials, with contours from -0.025 to 0.025 ms“^ with an interval of 0.005 ms“^, 
(a), and -0.005 to 0.005 ms“^ with an interval of 0.0005 ms“^, (b). 


To further evaluate the performance of the four model runs, the Root Mean 
Squared Error (RMSE) was determined for each run in relation to the analytic solution 
(see Table VI) for the four variables: vr, u, w, and 6. Once again, with the domain 
of interest defined and the numerical and analytic solutions on matching grids a 
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Table V. Case 3: Linear Non-Hydrostatic Mountain. Comparison of modified and 
unmodified CGI and CG2 using 360 m (in x) and 300 m (in z) resolution after 5 h. 


Model 

CGI x-z 

CGI x-dx 

CG2 x-z 

CG2 x-(T2 

max 

0.2281x10-^ 

0.2059x10-6 

0.2276x10-6 

0.2045x10-6 

tt ' ■ 

mm 

-0.2631x10-6 

-0.2628x10-6 

-0.2599x10-6 

-0.2654x10-6 

Umax (mS“^) 

0.8257x10-2 

0.8194x10-2 

0.8229x10-2 

0.8191x10-2 

Umin (mS"^) 

-0.7163x10-2 

-0.6910x10-2 

-0.7120x10-2 

-0.6903x10-2 

Wmax (mS"^) 

0.6035x10-2 

0.6696x10-2 

0.6035x10-2 

0.6694x10-2 

Wmin (mS“^) 

-0.6037x10-2 

-0.7989x10-2 

-0.6037x10-2 

-0.7988x10-2 

O'max (K) 

0.2962x10-2 

0.4275x10-2 

0.2969x10-2 

0.4262x10-2 

OLn (K) 

-0.2846x10-2 

-0.3987x10-2 

-0.2846x10-2 

-0.3991x10-2 

CPU time (s) 

0.6166x10^ 

0.9376x10^ 

0.1720x106 

0.1980x106 


bootstrap (random sampling) method, using 10,000 iterations, was used to calculate 
the 95 % Cl [11]. The resulting conhdence intervals again indicate that the differences 
between the state variable RMS errors are signihcant, since the ranges do not overlap. 
Both CGI x-z and CG2 x-z exhibited a smaller error for all four variables of interest, 
with the most signihcant difference in RMSE values for w and 9. Similar to the 
previous case, taking a closer look at the w RMSE with respect to Wmax shows that 
CGI x-z has a 2.17% error as compared to 20.25% for CGI x-a^ and that CG2 x-z 
has a 2.17% error as compared to 20.26% for CG2 x-a^, which is even larger than 
observed in the linear hydrostatic case. 

Similar to the previous case, in order to verify that a steady-state solution for 
the linear non-hydrostatic mountain case was achieved, the momentum hux was again 
derived using Eq. (5.1). [12] For this case, the analytic non-hydrostatic momentum 
hux is given by: 


= —0.457m^(;2) 


where denotes the analytic non-hydrostatic momentum hux and denotes 
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the analytic hydrostatic momentum flux (Eq. (5.2)). [13] The momentum flux was 
then normalized by m{z)/m^^ {z), and will be the value discussed in this thesis. In 
Figure 14, the normalized momentum flux for all four model runs, using 360 m (in x) 
and 300 m (in z) resolution, was plotted for 1 h, 2 h, 3 h, 4 h, and 5 h. The resulting 
momentum flux for this case resembled the same pattern as seen in case 2. Figure 14 
indicates that all four models converge to a steady state solution after 5 h. Of the 
four model simulations, the CGI x-z and CG2 x-z models yielded results that were 
far better than the x-a^ models, with values between 0.95 and 1.0 versus between 
0.85 and 1.18 as seen for the x-a^ runs. 

2. Conclusions 

All four model runs for the linear non-hydrostatic mountain case did develop 
a steady-state mountain wave over a single peak indicating that the model dynamics 
are resolving the scenario. Additionally, the models converged to nearly identical 
patterns, but with varying oscillation intensities. Similar to the results seen in case 2, 
for both set 1 and set 2 of the Navier-Stokes equations, the x-cr^ coordinates performed 
worse than their x-z counterparts (reflected in the RMSE, the normalized momentum 
flux , and most pronounced in w with an approximately nine times larger RMSE 
with respect to Wmax ) and did show a noticeable increase in computational expense. 
These results further indicate that when the model dynamics are resolved using purely 
explicit time integration methods that there is no significant benefit in using x-cTz over 
x-z and that there is a significant degradation, but the results are promising since the 
models are converging to a fairly representative steady-state solution. Although the 
Alters and boundary conditions are designed to be independent of z and az, it cannot 
be ruled out that a special treatment could be explored. With more research into the 
Alters and boundary conditions, the degradation could be made minimal enough that 
x-az is worth using with semi-implicit methods. 


63 



Table VI. Case 3: Linear Non-Hydrostatic Mountain RMSE. Root-mean-squared er¬ 
rors for the four primary state variables, for both the modihed and unmodihed codes, 
after 5 h using 360 m ( in x) and 300 m (in z) resolution and 10th order polynomials. 


Model 

71 

95% Conhdence Interval 

CGI x-z 

1.5939x10-^ 

+ 1.0987x10-4^^ 

- 1.1296x10-40 

CGI X-(Tz 

2.2354x10-® 

+ 1.8145x10-40 

- 1.8058x10-40 

CG2 x-z 

1.6391x10-® 

+ 1.1292x10-40 

- 1.1720x10-40 

CG2 x-(Tz 

2.2931x10-® 

+ 1.8913x10-40 

- 1.8253x10-40 

Model 

u (ms-^) 

95% Conhdence Interval 

CGI x-z 

4.6827x10-^ 

+ 3.5124x10-0 

- 3.5127x10-0 

CGI x-(T2 

5.3962x10-^ 

+ 4.3611x10-0 

- 4.3393x10-0 

CG2 x-z 

4.8941x10-^ 

+ 3.6237x10-0 

- 3.6915x10-0 

CG2 x-a. 

5.3838x10-^ 

+ 4.3711x10-0 

- 4.3818x10-0 

Model 

w (ms-^) 

95% Conhdence Interval 

CGI x-z 

1.3125x10-4 

+ 1.2912x10-0 

- 1.3101x10-0 

CGI x-a^ 

1.6180x10-® 

+ 1.2733x10-0 

- 1.2822x10-0 

CG2 x-z 

1.3105x10-4 

+ 1.3205x10-0 

- 1.3154x10-0 

CG2 x-a^ 

1.6180x10-® 

+ 1.2727x10-0 

- 1.2691x10-0 

Model 

0(K) 

95% Conhdence Interval 

CGI x-z 

1.6683x10-4 

+ 1.3953x10-0 

- 1.4053x10-0 

CGI x-a^ 

5.3367x10-4 

+ 4.7779x10-0 

- 4.7821x10-0 

CG2 x-z 

1.7285x10-4 

+ 1.4982x10-0 

- 1.4429x10-0 

CG2 x-cr^ 

5.3364x10-4 

+ 4.7514x10-0 

- 4.7686x10-0 
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ii) 




iii) 


iv) 


Figure 14. Case 3: Linear Non-Hydrostatic Mountain. Normalized momentum flux 
for the resolution of 360 m (in x) and 300 m (in z), at times 1 h, 2 h, 3 h, 4 h, and 
5 h for the four model runs: (i) CGI x-z, (ii) CGI x-a^, (iii) GG2 x-z, and (iv) GG2 
x-cr^. 
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VI. CONCLUSIONS AND 
RECOMMENDATIONS 

The results from the three test cases yielded promising results. The hrst of 
which was that the coordinates functioned properly and did in fact reduce to x-z 
coordinates in the absence of terrain, as demonstrated in case 1. With the introduc¬ 
tion of terrain, all four models converged to nearly identical steady state patterns for 
both case 2 and case 3 respectively, but had varying oscillation intensities between 
coordinate systems. Additionally, both case 2 and case 3 x-a^ model runs did de¬ 
velop a steady-state mountain wave over a single peak (reflected in the normalized 
momentum flux) indicating that the model dynamics are resolving the scenarios for 
a linear hydrostatic mountain and a linear non-hydrostatic mountain. In general, the 
CGI x-a^ and CG2 x-a^ models performed worse than their x-z counterparts. The 
x-az models had higher RMS errors (as large as nine times greater RMSE values with 
respect to associated maximum vertical velocities), which were observed predomi¬ 
nantly in intensity values and not in placement of steady state features. All three 
cases for x-a^ also showed a noticeable increase in computational expense, due to 
the additional calculations and variables required by the coordinate transformation. 
These results indicate that though x-az coordinates are not as accurate or efficient 
as x-z and that there is a signihcant degradation, with further hne-tuning of the 
model environment these issues could be made minimal enough that it justihes their 
use in semi-implicit methods, especially in the vertical, as is already done by various 
operational mesoscales models. 

Past research has demonstrated the utility of using Continuous Galerkin meth¬ 
ods, in a x-z framework, for resolving computational fluid dynamics using both fully 
explicit and semi-implicit time integration methods. This thesis brought Giraldo’s 
2-D (x-z slice) mesoscale Non-Hydrostatic model one step closer to evaluating and 
exploiting the full strength of x-az coordinates by transforming the Navier-Stokes 
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Equations and testing their impacts using fully explicit time integration. Now that 
the x-(Tz models are functional, the next stage of development will be to implement a 
semi-implicit time integration method in the vertical (1-D) and compare the results 
against the 2-D semi-implicit x-z model. Without the work done in this project to 
transform the x-z eqnations to x-a^, one could not construct semi-implicit methods in 
the vertical since in x-z, the terrain is conpled to box coordinates and so the coordi¬ 
nates cannot be decoupled. Additionally, the mathematically machinery outlined in 
this thesis can be used to transform any eqnation set to any other coordinate system. 

The valne of this study is far reaching in determining the usefulness of ap¬ 
plying a specihc coordinate system in the fnture when developing meteorological and 
oceanographic models for the U.S. Naval Research Laboratory (NRL) by constituents 
at the Naval Postgradnate School. In addition, the snccessful conversion of the non¬ 
hydrostatic x-z models to x-a^ will allow for the straightforward extension of these 
models to global non-hydrostatic models, since a will represent the height of the 
model and x will then represent the spherical shell at each valne of a. 



APPENDIX. COEFFICIENTS FOR RK35 

METHOD 


CoefRcient Value 


stage 1: 

Uq 

1 


0.377268915331368 

stage 2: 


0 

Ul 

1 

Po 

0.377268915331368 

stage 3: 

tto 

0.355909775063327 

CXi 

0 

Oi2 

0.644090224936674 

Po 

0.242995220537396 

stage 4: 

Uq 

0.367933791638137 

CXl 

0 

Oi2 

0 

Us 

0.632066208361863 

Po 

0.238458932846290 

stage 5: 

ao 

0 

tti 

0 

Oi2 

0.237593836598569 

«3 

0 


0.762406163401431 

Po 

0.287632146308408 
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